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Abstract 

Various artificial compressibility methods for calculating the three-dimensional 
incompressible Navier-Stokes equations are compared. Each method is de- 
scribed and numerical solutions to test problems are conducted. A compari- 
son based on convergence behavior, accuracy, and robustness is given. 


1 Introduction 

The difficulty in computing solutions to the incompressible Navier-Stokes sys- 
tem of PDEs lies in satisfying the divergence-free velocity condition. Artificial 
compressibility methods, developed by A. Chorin [1] , provide a mechanism 
to march in pseudo-time towards the divergence-free velocity field such that 
mass and momentum are conserved in the pseudo steady-state. 

The classical artificial compressibility method transforms the mixed ellip- 
tic/parabolic type equations into a system of hyperbolic or parabolic equa- 
tions in pseudo-time, which can be numerically integrated. The method has 
been generalized to curvilinear coordinates and used for various applications, 
Kiris et. al. [3], 

Since the publication of Chorin's original paper many alternative forms 
of artificial compressibility have been developed. These methods include a 
generalized preconditioning matrix to equalize the wave speeds and the use of 
differential preconditioning, Turkel and Radespiel [7], as well as the addition 
of artificial viscosity such as an artificial Laplacian of pressure term in the 
continuity equation, Shen [6]. We present a direct comparison of four different 
versions of the artificial compressibility method on a series of test problems. 


2 Incompressible Equations and 
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The governing equations for incompressible, constant density and constant 
viscosity flow written in conservative form in generalized coordinates with 
the density absorbed by the non-dimensionalization of the pressure term are 
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Classical and Generalized Artificial Compressibility 


( 1 ) 


The classical artificial compressibility formulation is derived by introduc- 
ing an artificial compressibility relation in the continuity equation. This is 
achieved by adding a preconditioned pseudo-time derivative of the primitive 
variables Q to equation 1. The generalization of this approach is to begin 
with the conserved variables W = (p, pu, pv, pw) T and use the chain rule to 
derive the generalized preconditioned pseudo-time derivative. The classical 
preconditioning matrix, T c , and the generalized preconditioning matrix, F p , 
are 
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where /? > 0 is the artificial compressibility parameter. 

Artificial Dissipation 

To assist in the dissipation of spurious pressure waves we add an artificial 
Laplacian of pressure term to the right-hand side of the classical artificial 
compressibility continuity equation. This term is scaled by a parameter e 
and has the affect of adding a second difference artificial dissipation term to 
the continuity equation. 


Differential Preconditioning 


The artificial dissipation described above manipulates the physical dissipation 
properties of the PDE system. Alternatively we can manipulate the convective 
properties of the system. Following Turkel and P.adespiel [7] a pseudo-time 
derivative of the Laplacian of pressure can be added to the continuity equation 
along with the standard pseudo-time derivative of pressure. This term is 
also scaled by e and will have an effect of propagating the low-frequency 
components of error more quickly than the high-frequency components which 
will be dissipated by the discretization scheme. 
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3 Numerical Results 

The INS3D code [4], [5], [3] has been adapted to include each of the artificial 
compressibility methods described. An implicit line symmetric Gauss-Seidel 
relaxation scheme is used with fully-implicit boundary conditions. Iterations 
are performed until the residual of the no nlin ear system has been reduced nine 
orders of magnitude in the l 2 norm. Results for 0 = 1, 10, 100 and CFL = 1 
and CFL — 1000 sire provided. For the artificial dissipation method e values 
of 1.0 -2 , 1.0 -1 , 1.0 +0 me used to scale the Laplacian term. The differential 
preconditioned method uses values of e = 1.0 -1 , 1.0 +o , 1.0 +1 . For each test- 
case the inlet velocity is specified and a constant pressure is enforced at 
the outlet. P = 1,2, 3, 4 denotes the version of the artificial compressibility 
method. 

Test 1: Inviscid Flow in a Square Duct 

Each method is used to calculate the inviscid flow in a square duct with 
dimensions 10 x 1 x 1 non-dimensional units. The exact solution is Q = 
(0, 1,0, 0) T . A grid of dimension 33x9x9 is used. Table 1 displays the number 
of iterations required. Each of the methods computed the correct solution up 
to double precision. For 0 = 1 the generalized preconditioned method has the 
best convergence. For (3 > 1 , the classical has the best convergence rate with 
the exception of CFL = 1000 where the differential preconditioning scheme 
is slightly better. 


Table 1 . Inviscid square duct: Number of iterations for residual reduction of 9 
orders of magnitude in the discrete L 2 norm. 
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Test 2: Viscous Flow in a Circular Pipe 

A simple viscous flow in a circular pipe of radius one and length ten is com- 
puted. A Reynolds number of 1000 is used for which an exact solution is 
d^riv^d. Grids of dimsiisioii IT x 9 x 9, 33 x IT x IT, snd 65 x 33 x 33 sx q nssd. 
Each method is verified to produce second order accurate results for 0 = 10. 
Figure 1 plots the normalized l 2 residual for varying 0 and CFL = 1000. 
Table 2 displays the number of iterations required to converge on the finest 
mesh. The third and fourth methods fail to converge for low CFL numbers 
and the artificial dissipation scheme is especially sensitive to the e parameter. 
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For high CFL and high /3 the differential preconditioning becomes effective, 
but a moderate f3 must be used for accuracy purposes. 




Fig. 1. Viscous Straight Pipe (Grid 65 x 33 x 33, Re = 1000, CFL = 1000): 
Comparison of the convergence between preconditioning methods for (3 = l(left) 
and /3 = 10 (right), o P = 1; □ P = 2; x P = 3, 6 = 1.0~ 2 ; + P = 3,€=1.0“ 1 ; 
* P = 3,e= 1.0+°; o P = 4,e= 1.0 -1 ; A P = 4, e = 1.0 +o ; * P = 4,e = 1.0 +1 . 


Table 2. Viscous Straight Pipe: Number of iterations for residual reduction of 9 
orders of magnitude in the discrete L 2 norm. 
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Test 3: Viscous Flow in an L-shaped Duct 

A more complicated viscous flow in a square duct with a 90° bend is used for 
the final test. The geometry used is described in Humphrey [2], where exper- 
imental results were obtained for Reynolds number 790. A grid of dimension 
65 x 33 x 33 is used. Figure 2 plot the residual for varying (3 and CFL = 1000. 
Table 3 displays the number of iterations required. The symbol * * ** denotes 
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the method failed to converge.Figure 3 displays a comparison of the different 
computed solutions with the experimental data for ft = 1 , 10, 100 at 9 = 90° 
plane of the curved duct. Robustness is an issue for the third and fourth 
schemes. Comparing the experimental data with the computed solutions we 
observe that using ft > 10 leads to poor solution accuracy and ft — 1 is the 
most accurate. 




Fig. 2. Viscous Square Duct with 90° bend (Grid 65 x 33 x 33, Re = 790, CFL = 
1000): Comparison of the convergence between preconditioning methods for ft = 1 
(left) and ft — 10 (right), o jP = 1; □ P — 2; x P = 3,e = 1.0~ 2 ; + P — 
3, e = 1.0 1 ; * P = 3,e = 1.0 +o ; o P = 4,e = ID" 1 ; A P = 4,e = 1.0+°; 
* P = 4,e= 1.0+ 1 . 


4 Conclusion 

Four variations of the artificial compressibility method have been imple- 
mented and compared on a series of test problems. The classical and gen- 
eralized artificial compressibility methods have almost identical convergence 
rates on each test case for every combination of the CFL and ft parameters. 
The artificial dissipation and differential preconditioning methods were not 
able to converge for all parameters on certain problems showing a lack of 
robustness for these methods. High values of ft lead to poor accuracy for all 
the methods considered. For moderate values of 1 < ft < 10 the classical and 
generalized versions appear to be the most accurate. These two versions will 
be evaluated for more complicated engineering applications. 
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Table 3. Viscous Square Duct with 90° bend: Number of iterations for residual 
reduction of 9 orders of magnitude in the discrete L 2 norm. 
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Fig. 3. Comparison of the solution at 9 = 90° and p = 1 (left), (3 = 10 (center), 
and P = 100 (right), o Exper.\ □ P = 1; A P = 2; x P = 3; • P = 4. 
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Abstract 


Various artificial compressibility methods for 
calculating the three-dimensional incompress- 
ible Navier-Stokes equations are compared. Each 
method is described and numerical solutions 
to test problems are conducted. A comparison 
based on convergence behavior, accuracy, and 
robustness is given. 
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Introduction 


The difficulty in computing solutions to the 
incompressible Navier-Stokes system of .PDEs 
lies in satisfying the divergence-free velocity 
condition. Artificial compressibility methods, 
developed by A. Chorin [1], provide a mech- 
anism to march in pseudo-time towards the 
divergence-free velocity field such that mass 
and momentum are conserved in the pseudo 
steady-state. 

The classical artificial compressibility method 
transforms the mixed elliptic/parabolic type equa- 
tions into a system of hyperbolic or parabolic 
equations in pseudo-time, which can be nu- 
merically integrated. The method has been 
generalized to curvilinear coordinates and used 
for various applications, Kiris et. al. [2]. 
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Since the publication of Chorin's original paper 
many alternative forms of artificial compress- 
ibility have been developed. These methods in- 
clude a generalized preconditioning matrix to 
equalize the wave speeds and the use of dif- 
ferential preconditioning, Turkel and Radespiel 
[3], as well as the addition of artificial viscos- 
ity such as an artificial Laplacian of pressure 
term in the continuity equation, Shen [4]. We 
present a direct comparison of four different 
versions of the artificial compressibility method 
on a series of test problems. These problems 
include inviscid flow in a square duct, viscous 

A - 

flow in a circular pipe, and viscous fiow in a 
square duct with strong curvature. 
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Governing Equations 


The governing equations for incompressible, 
constant density and constant viscosity flow 
written in conservative form in generalized co- 
ordinates with the density absorbed by the nond 
mensionalization of the pressure term are 
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Artificial Compressiblity Methods 


The artificial compressibility formulation is de- 
rived by introducing an artificial compressibility 
relation, 



( 2 ) 


where p is the pressure and (3 > 0 is the artificial 
sound speed. Using this relation we may add 
a preconditioned pseudo-time derivative of the 
primitive variables Q to equation 1. Four forms 
of artificial compressiblity will be discussed. 


• Classical Artificial Compressiblity 


• Generalized Artificial Compressibility 


• Artificial Dissipation 


• Differential Preconditioning 
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Classical Artificial Compressiblity 


The classical artificial compressibility method 
uses equation 2 only in the continuity equation. 


This leads to the following preconditioned sys- 
tem of equations, 


r dQ 8Q , d{E t - EV) 

fc ^ + Im ^ + 5T 




The classical preconditioning matrix is given 
by, 
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Generalized Artificial Compressiblity 


To generalize the above approach we start with 
the conserved variables W = (p, pu, pv, pw) T . 
Using the chain rule and equation 2 we obtain, 

dW _ dWdQ _ n dQ 
~dQ^¥r~ P lh 1, 

where 

A 0 0 0 
110 0 
| 0 1 0 
f 0 0 1 

Substituting r p for V c we obtain the general- 
ized precondition system of equations, 
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Artificial Dissipation 


Introducing a finite sound speed into the in- 
compressible equations creates artificial pres- 
sure waves which must be propagated out of 
the solution domain in pseudo-time. 


Alternatively the artificial pressure waves can 
be dissipated by adding an artificial Laplacian 
of pressure to the modified continuity equation 

in equation 3. Where the viscous fluxes Ej' are 
replaced by, 
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Differential Preconditioning 


The classical artificial compressibility method 
uses a standard pseudo-time derivative to con- 
vert the artificial pressure waves out of the 
domain. So both high and low frequency er- 
rors are convected at constant wave speeds. 
Alternatively, the time derivative of the Lapla- 
cian of pressure term can be added to the first 
equation in the system. This forces the low 
frequecy errors to travel at faster speeds then 
there high frequency counterparts. The high 
frequency errors will be damped by the dis- 
cretization and relaxation schemes. 


This produces a system of the form, 
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Numerical Examples 

The INS3D code has been adapted to include each of the artifi- 
cial compressibility methods described. An implicit line symmetric 
Gauss-Seidel relaxation scheme is used with fully-implicit boundary 
conditions. Iterations are performed until the residual has been 
reduced nine orders of magnitude in the l 2 norm. Results for 
f 3 = 1,10,100 and CFL = 1 and CFL = 1000 are provided. For 
the artificial dissipation method e values of 1.0 -2 and 1.0 -1 are 
used to scale the Laplacian term. The differential preconditioned 
method uses values of e = 1.0 -1 and 1.0 +1 . For each test-case the 
inlet velocity is specified and a constant pressure is enforced at the 
outlet. 

• Inviscid Square Duct 
e viscous Circular Pipe 

• Viscous Square Duct with Strong Curva- 
ture 
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Inviscid Square Duct 


Each method is used to calculate the inviscid 
flow in a square duct with dimensions 10 x 1 x 1 
non-dimensional units. The exact solution is 
Q — (0, 1, 0, 0) r . A grid of dimension 33 x 9 x 9 
is used. The table below displays the num- 
ber of iterations required. Each of the meth- 
ods computed the correct solution up to dou- 
ble precision. For (3 = 1 the generalized pre- 
conditioned method has the best convergence. 
For p > 1, the classical has the best conver- 
gence rate with the exception of CFL = 1000 
where the differential preconditioner scheme is 
slightly better. P — 1,2, 3, 4 denotes the ver- 
sion of the artificial compressibility method. 
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Grid 33 x 9 x 9, Re — 1000, CFL = 1000: Plot 
of normalized l 2 norm residual (3= 10 


15 





DOO, CFL = 1000: Plot 
residual f3 = 100 
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Viscous Circular Pipe 


Viscous flow in a circular pipe of radius one 
and length ten is computed. A Reynolds num- 
ber of 1000 is used for which an exact solution 
is derived. A Grid of dimension 65 x 33 x 33 
is used. The table below displays the number 
of iterations required to converge on the finest 
mesh. The third and fourth methods fail to 
converge for low CFL numbers and the artifi- 
cial dissipation scheme is especially sensitive to 
the e parameter. For high CFL and high /? the 
differential preconditioning becomes effective, 
but a moderate (3 must be used for accuracy 
purposes. 
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Iterations 

Grid 65 x 33 x 33, Re = 1000, CFL = 1000: 
Plot of normalized l 2 norm residual (3=1 
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Log 10 (||R n ||/||R°||) 



Grid 65 x 33 x 33, Re = 1000, CFL = 1000: 
Plot of normalized l 2 norm residua! f3 = 10 
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Grid 65 x 33 x 33, Re = 1000, CFL = 1000: 
Plot of normalized l 2 norm residual /? = 100 
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Viscous Square Duct with Strong 

Curvature 


Viscous flow in a square duct with a 90° bend 
is used for the final test. The geometry used is 
described in Humphrey [5], where experimen- 
tal results were obtained for Reynolds number 
790. A grid of dimension 65 x 33 x 33 is used. 
The table below displays the number of itera- 
tions required to converge. The symbol * * ** 
denotes the method failed to converge. Ro- 
bustness is an issue for the third and fourth 
schemes. Comparing the experimental data 
with the computed solutions we observe that 
using (3 > 10 leads to poor solution accuracy 
and ( 3=1 is the most accurate. 
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Grid 65 x 33 x 33, Re = 1000, CFL = 1000: 
Plot of normalized l 2 norm residual (3 = 1 
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Grid 65 x 33 x 33, Re = 1000, CFL = 1000: 
Plot of normalized l 2 norm residual (3 = 10 
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Grid 65 x 33 x 33, Re = 1000, CFL = 1000: 
Plot of normalized l 2 norm residua! (3 = 100 
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Conculsion 


Four variations of the artificial compressibility 
method have been implemented and compared 
on a series of test problems. The classical 
and generalized artificial compressibility meth- 
ods have almost identical convergence rates 
on each test case for every combination of the 
CFL and ,3 parameters. The artificial dissipa- 
tion and differential preconditioning methods 
were not able to converge for all parameters on 
certain problems showing a lack of robustness 
for these methods. High values of /3 lead to 
poor accuracy for all the methods considered. 
For moderate values of 1 < f3 < 10 the classical 
and generalized versions appear to be the most 
accurate. These two versions will be evaluated 
for more complicated engineering applications. 
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